############################################################################
############################################################################
####Non-state Atrocities in Capital Cities - ZINB Global Sample Analysis####
############################################################################
############################################################################
library(MASS) 
library(pscl)
library(foreign)
library(stargazer)
library(mvtnorm)
library(plyr)

# Set working library
setwd("~/Data/Global Analysis/")
# Read in main data
main.data <- read.dta("full.grid.upd.dta")

##############
###Baseline###
##############
###Model 1: baseline
at.zinb.1 <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull
                      + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                        year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                      |lagcivconflagtemp + loglagcellarea + loglagpop + loglagttime, 
                      dist = "negbin", data = main.data)
summary(at.zinb.1)
AIC(at.zinb.1)

#############################
###Add Inflation Variables###
#############################
at.zinb.2 <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull
                      + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                        year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                      |lagcivconflagtemp + lagincidentsumfull
                      + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                      dist = "negbin", data = main.data)
summary(at.zinb.2)
AIC(at.zinb.2)

######################
###Add Urbanization###
######################
at.zinb.3 <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban
                      + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                        year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                      |lagcivconflagtemp + lagincidentsumfull + lagurban
                      + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                      dist = "negbin", data = main.data)
summary(at.zinb.3)
AIC(at.zinb.3)

########################
###Add Oil Production###
########################
at.zinb.4 <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                      + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                        year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                      |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                      + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                      dist = "negbin", data = main.data)
summary(at.zinb.4)
AIC(at.zinb.4)

###Export to latex
stargazer(at.zinb.1, at.zinb.2, at.zinb.3, at.zinb.4)
stargazer(at.zinb.1, at.zinb.2, at.zinb.3, at.zinb.4, zero.component=TRUE)
AIC(at.zinb.1, at.zinb.2, at.zinb.3, at.zinb.4)


##################################################
####Vuong Tests -- Negative Binomial Comparison###
##################################################
##Models 1 and 2
nb1 <- glm.nb(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull
              + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08, 
              data = main.data)
#Compare NB and ZINB Models 1 and 2
vuong(nb1, at.zinb.1)
vuong(nb1, at.zinb.2)

##Model 3
nb3 <- glm.nb(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban
              + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop +  loglagttime + loglagcellarea + year96 + year97+
                year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08, 
              link = log, data = main.data)
#Compare NB and ZINB Model 3
vuong(nb3, at.zinb.3)

##Model 4
nb4 <- glm.nb(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
              + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop +  loglagttime + loglagcellarea + year96 + year97+
                year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08, 
              link = log, data = main.data)
#Compare NB and ZINB Model 4
vuong(nb4, at.zinb.4)
